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We explore the superfluid phases of a two-component Fermi mixture with hybridized orbitals in 
optical lattices. We show that there exists a general mapping of this system to the Lieb lattice. 

By using simple multiband models with hopping between s and p-orbital states, we show that 
superfluid order parameters can have a 7r-phase difference between lattice sites, which is distinct 
from the case with hopping between s-orbitals. If the population imbalance between the two spin 
species is tuned, the superfluid phase may evolve through various phases due to the interplay between 
hopping, interactions and imbalance. We show that the rich behavior is observable in experimentally 
realizable systems. 
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I. INTRODUCTION 

Multiband effects are important in understanding a va¬ 
riety of quantum many-body phenomena such as high- 
temperature superconductivity, fractional quantum Hall 
phases, and topological matter in general mM- In the 
context of ultracold gases, excellent opportunities of ex¬ 
perimentally exploring multiband phenomena in a con¬ 
trolled way have emerged in the recent years. Optical 
lattices naturally have multiple bands, and quantum sys¬ 
tems in the excited bands of optical lattices have been 
experimentally realized. Muller et al. [5] transferred ul¬ 
tracold bosons into the p-band of the lattice and observed 
how coherence was established between atoms. Recently 
Zhai et al. prepared bosons in the d-band [6]. Closely 
related to the work presented in this article are the ex¬ 
periments I3i] exploring excited band condensates in 
a bipartite optical lattice, and m, studying pairing be¬ 
tween different parity orbital fermions. Bosons in the flat 
band of a Lieb lattice were realized very recently m- 

These experimental possibilities have insprired consid¬ 
erable amount of theoretical work on multiband effects in 
optical lattices. Naturally, in other contexts the volume 
of work on multiband effects is much larger; here we men¬ 
tion only examples of results related to ultracold gases. 
For example, it has been argued that since higher bands 
of an optical lattice typically have larger bandwidths, 
higher critical temperatures for anti-ferromagnetic order¬ 
ing may be realized m- The idea is probably mentioned 
in many places and is based on the fact that in perturba¬ 
tion theory the Hamiltonian ends up having a prefactor 
tunneling squared over coupling. In the p-band tunneling 
can be larger so the chracteristic energy over temperature 
scale can be higher in absolute sense. (Quote from Wu 
et al. [12] : “We also show that in the strongly correlated 
regime the Neel temperature for p band antiferromag¬ 
netism is 2 to 3 orders of magnitudes higher than that of 
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s band, which is much more promising to be attained in 
cold atom experiments.”) It has also been suggested by 
Dutta et al. na that strong interactions and multiband 
effects can give rise to self-assembly of non-trivial lattices 
for topological insulators. Topological semimetals and 
chiral superfluidity with s-wave interactions have been 
predicted nana in multi-orbital models where orbitals 
with different symmetries interact. There is also exten¬ 
sive theory literature on both p-band bosons [MIS] and 
fermions [20| with many studies focusing on the strong 
coupling regime 1211 [22]. Furthermore, oscillating order 
parameters in fermionic systems have been predicted due 
to coupling between s- and p orbitals [23] , or in pure p- 
orbital systems [24l [25] . 

Motivated by these advances, in this article we explore 
the physics of attractively interacting two-component 
fermions with multiple bands. In particular, we wish to 
understand how an unequal number of different fermionic 
species and the tunneling properties of p-orbitals influ¬ 
ence the formation of s-wave pairing order parameters in 
such systems. 

We first solve a simplified model with just two 
sites which indicates possibilities of different superfluid 
phases, some with a spatially varying order parameter 
phase factor. We then demonstrate that in many re¬ 
spects, the results from this toy model are realized by a 
bipartite lattice where s and p-orbitals in different sub¬ 
lattices hybridize [26]. Such a lattice has been experi¬ 
mentally demonstrated by Wirth et al. [ 7 ] who studied 
Bose-Einstein condensation on the excited bands of such 
a lattice and found non-trivial ordering of the conden¬ 
sate phase. This ordering is due to the fairly complex in¬ 
terplay between tunneling properties of different orbitals 
and on-site interactions between atoms. 

We outline the expected phase diagram for fermions 
at the mean-field level and find a possibility of a tt- 
phase superfluidity where the sign of the order parameter 
varies between sublattices. Such possibility was raised by 
Iskin m in the context of a checkerboard lattice, but it 
turned out that in that system this possibility was not re¬ 
alized. Somewhat related phenomena have also been dis¬ 
cussed in studies exploring EELO phases in lattices [28L 
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SB, in multiorbital models [23l[25], or in two-dimensional 
systems without a lattice [32] . 

This paper is organized as follows: We start, in Sec. [HI 
by discussing the lattice we consider in some detail and 
we especially elaborate on the sign changes that occur 
in the hopping parameters, due to the different orbital 
states. Consequently, in Sec. |ml we introduce interac¬ 
tions and the possibility of a pairing instability. We first 
study a simplified, dispersionless model to examine what 
kind of pairing instabilities can occur. Subsequently, in 
Sec. IV, we study a more complex system, resembling an 
experimentally realizable system. The numerical results 
are presented in Sec. V. Finally in Sec. VI, we conclude 
with a summary and discussion. 


II. SIMPLE MODEL WITH HYBRIDIZED 
ORBITALS 



FIG. 1: The lattice we study in one dimension (a), where 
on the B lattice sites the p orbitals are sketched. Hopping 
within the specified unit cell is denoted by a full line, whereas 
hopping to the neighboring unit cell is denoted by a dashed 
line. In (b) the mapping from the s-p lattice to the lattice 
with s orbitals only is shown. 


We consider fermions in two different (pseudo) spin 
states, t and occupying bipartite lattices with sublat¬ 
tices A and B. Here, the fermions on sites of the A sub¬ 
lattice are in the 5-orbital state, whereas the particles in 
the B sublattice occupy a p-wave orbital state, where the 
number of the different p-orbitals is equal to the dimen¬ 
sion of the lattice. In this section, we study first non¬ 
interacting fermions occupying a one-dimensional and 
consequently a two-dimensional lattice, where we denote 
the orbitals on the B sublattice by p or p^ and Py^ re¬ 
spectively. Our goal here is to elucidate how systems 
with hybridized s and p orbitals are connected to and 
differ from systems with s orbitals only. 


A. One dimension 


Due to the different orbital states on the two sublat¬ 
tices, the hopping coefficients are also different for parti¬ 
cles moving in opposite directions. To be more specific, 
we first focus on the one-dimensional case. There, the 
hopping coefficient for a particle moving from an A to 
a neighbouring B site in the positive x-direction, 
has an opposite sign to the one in the opposite direc¬ 
tion, which is due to the odd parity of the p orbital. 
The hopping coefficients for the particles moving back 
from the B to the A sublattice are the same, 
and , since they correspond to the same overlap 

integral. Apart from the sign difference the hopping coef¬ 
ficients are the same, = t. Thus, the nearest 

neighbour hopping Hamiltonian for the one-dimensional 
case reads 




CinV’a.n “ ^a,n+lVa,n + h-C. 


lP\ 


( 1 ) 


where creates a (pseudo) spin state \a) fermion with 
orbital j at unit cell n and a unit cell contains one A and 
one B site, see FigQa). The first term thus describes 


hopping within a unit cell, while the second describes 
hopping across unit cells. 

It is instructive to Fourier transform the above Hamil¬ 
tonian Eq. 0, in a few steps 

H]F = -ty [+ h.c.' 

k,cr 

k,cr 

= -t y]2i sinikd) > (2) 

k,a 

where d is the lattice spacing, which is taken equal to one 
here. A transformation in the fermionic operators was 
made in the last line, '0^/. = exp[—^ and = 

k,k- 


B. Mapping p to s orbitals 

Although the Hamiltonian in Eq. 0 and its Eourier 
transform in Eq.([^ look quite different from the Hamil¬ 
tonian describing fermionic particles in a ID lattice with 
only 5-orbital sites, it turns out that these two are more 
similar than they seem. We show this connection by 
transforming the Hamiltonian with alternating hoppings 
to a hopping Hamiltonian without sign changes, see 
EigQb). Eirst, Eq.Q can be rewritten by splitting the 
summation over the unit cells in a sum over the even and 
odd unit cells, after which the Hamiltonian reads 

H]F = + ^2m+l,<T^L+l.<T 

m,a 

- i’tL+lAim+2,a + h-C. 
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where the summation over m runs over half of the values 
that n in Eq. 0 runs over. Consequently, the following 
unitary transformation can be used 

'>p2m,a = ^2m,a 

i’2m+l,a = -i’2m+l,a^ (3) 

with j = s,p. This yields for the Hamiltonian 

= -tYl [^2ln,a^2m,a + ^2m+l,a^fm+l.a 

m, a 

+^2m,<T^lm+l,<T + ^2TO+l,<T^lm+2,(T + ^.C. 

= -i E [^n!<T^n,a + + h-C. , 

n, cr 


which is the hopping Hamiltonian for particles in a ID 
lattice with only 5-orbital sites. In Fourier space the same 
transformation boils down to making a shift in the quasi 
momentum of the operators, ~ '^i- 7 r/( 2 d) cr* 
subsequently shifting all quasi momenta by +7r/(2(i) the 
Hamiltonian reads 

H¥’ = -i E [(1 - 

k,a 

+ (1 - 

= -t E2 cos{kd) > (4) 

k 



where the same transformation as before has been used 
in the last step. 

Although there thus exists a simple mapping between 
the s-p Hamiltonian and the 5 orbitals only lattice, from 
Eq.Q it is also clear that they are not exactly the same, 
the difference being the dispersions of the particles. The 
different dispersions result in a number of differences 
between the two systems, such as different momentum 
distributions and Fermi momenta. The latter can in 
turn give rise to different properties of possible super¬ 
fluid phases. 


C. Two dimensions 

In a two-dimensional system the lattice sites of the 
B sublattice contain two p-orbitals, Px and Py (denoted 
by X and y in sub- or superscripts, respectively), see 
FigJl a). The presence of two p-orbitals changes the 
hopping physics yet a bit more compared to the one¬ 
dimensional case. Namely, a particle that moves from 
the A to the B sublattice in the x direction, ends up 
in a Px orbital state. Since the overlap integral between 
the Px orbital and the neighboring 5-orbitals in the y- 
direction vanishes, this particle can only move along the 
x-direction. And similarly for particles moving from the 
A to the B sublattice in the y direction. In other words. 


FIG. 2: (a) The s-p lattice in two dimensions, where on the 
B lattice sites the px and py orbitals are sketched. Hopping 
within the specified unit cell is denoted by a full line and to 
the neighboring unit cells by a dashed line, (b) Illustration 
of the two decoupled lattices. On the left the full lattice is 
sketched and the two different hopping routes for the particles 
are denoted by full and dashed lines, which are shown as 
separate lattices in the middle and right, (c) mapping from 
the s-p lattice to the lattice with s orbitals only. 


^±x — ^±y — 0 parity of the p-orbitals. In 

the absence of interactions, this results in effectively two 
hopping sublattices, illustrated in Fig. [^b), which both 
have the Lieb lattice geometry. The hopping Hamilto¬ 
nian in two dimensions reads 


H 


2D 

K 


E 

n,cr;Q:G{cc,?/} 





cr,n-|-Q; ‘r cr ,n 


Vo 


h.c. 


(5) 


where a unit cell contains one site from the A and one 
from the B sublattice, see Figj^a), and the second term 
in the summation now describes hopping to the next unit 
cell both in the x and y direction, denoted by n + x and 
n y respectively. 

The Fourier transform of the hopping Hamiltonian for 
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the two-dimensional system is 


III. PAIRING 


= -t Y. 2i sin(fc„d) (ClfeC,;. - €]k^a,k) , 

/e,cr;aG{a:^,?/} 

( 6 ) 

where d is again the lattice spacing and the same trans¬ 
formation in the fermionic operators has been used as in 
the one-dimensional case Eq. <§• 


D. Mapping px and py to s orbitals 

Also in a two-dimensional system a mapping can be 
made from the s-p lattice discussed above to a lattice 
containing only s-orbital states. The mapping from a 
square lattice containing both s and p orbitals is to the 
Lieb lattice containing only s orbitals, see FiglUJc). Ac- 
tually, the two hopping sublattices both map to a Lieb 
lattice. Since they are uncoupled, the Hamiltonian in 
Eq.([^ can also be written as the sum of the two hopping 
Hamiltonians of the sublattices 

=Hki+Hk2, 


In the previous section we studied particles in bipartite 
lattices and looked in detail at the effect of the different 
orbitals on the hopping. To explore the role of parity 
further, in this section we include on-site interactions 
between the fermionic particles in the one-dimensional 
lattice. We specifically look at pairing instabilities that 
can arise due to attractive interactions and study the dif¬ 
ferent superfluid phases that occur in this system. This 
section is meant as an instructive example, since the no¬ 
tion of a superfiuid in one dimension is more complicated 
than the definition we use here. 

The interaction Hamiltonian reads 




ID 








( 7 ) 

where we consider attractive interactions, < 0, 

which in general can be different. We use a mean- 
field approximation and include Cooper pairs Aq = 

Hamiltonian reads 


where both Hki and Hk 2 have the same form as Eq.([^, 
but both now sum over half of the unit cells, see EigMl^. 
Just as before, the summation can be split in a sum over 
the even and odd unit cells, which for Hki reads 


= ijy -h hd 

= + (8) 

k 


Hki = -t 


E 

m,a;aE{x,y} 
A“t 


^a!2mC,2m 






'cr, 2 m+l'^a 


After using the same transformation for the fields as in 
the one-dimensional case, Eq.([^, the Hamiltonian reads 


Hki =-t Y 

m,a-,aE{x,y} 


^a,2m^a,2m + C!2m+lC,2m+l 


+C,2m+a<,2m 

= -^ E 


+ Cj2m+l)+aC,2m+l + h.C._ 
K\n'^a,n' + K]n'+a'^l,n' + ^-C. 


n',cr;aG{a:,y} 


where the summation over n' runs over half of the values 
that n runs over in the full 2D Hamiltonian, Eq. (|5|. The 
above Hamiltonian now describes particles hopping in 
a two-dimensional Lieb lattice with s-orbital sites only. 
The same transformation can be made for Hk 2 - As in the 
one-dimensional case, the dispersions describing particles 
in the lattice with both s and p orbitals are of sine form, 
see Eq. while particles in the Lieb lattice with only 
s orbitals have a cosine dispersion. The Lieb lattice has 
been studied in the context of cold gases as well [33l 
[34] , where interesting phenomena result from the lattice 
exhibiting flat dispersions, so-called flat bands. Recently, 
Bose-Einstein condensation of atoms in a Lieb lattice was 
experimentally studied by Taie et al. m- 


where is the chemical potential for fermions in spin 
state \a) and where the matrix Elves' in the Nambu basis 
with becomes 


^BCS = 



(9) 


with the average chemical potential /i = 
and where the possibility of having an imbalance in the 
population of the two spin components is included via a 
chemical potential difference h = (/i^ — /i4,)/2. Eor the 
spin-down sector we used S-k = sin(—/c) = —Sk for the 
dispersions in the hopping Hamiltonian Eq. © and we 
set d = 1. In order to write the total Hamiltonian using 
matrix multiplication, the spin-down fields have been in¬ 
terchanged. In the usual BCS theory this would result 
in extra terms Sk — ES], whereas here only the 2/i^ 
in Eq.(|^ stems from interchanging fermionic fields. This 
is due to the alternating signs for the hoppings, meaning 
that the two dispersions coming from interchanging spin- 
down fields cancel each other. It is assumed in the above 
Hamiltonian that there is no energy offset between the A 
and B lattice sites. 

To understand the different phases that can occur in 
this system, we first neglect the momentum dependen¬ 
cies of the particle dispersions by setting them all equal 
to one, 5/c = 1, for simplicity. By diagonalizing the above 
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Hamiltonian the four quasiparticle dispersions hwi are 
obtained, which apart from the usual |Aop and |Aip 
terms, now also contain mixed terms, such as AqAi. 
Consequently, the partition function can be obtained 
Z = Tr[exp(—ydi^)], from which in turn the thermody¬ 
namic potential can be calculated ft = — lnZ//3, where 
P = l/ksT is the inverse thermal energy, with ks Boltz¬ 
mann’s constant [35]. For our system the thermodynamic 
potential reads 




|AoP 

Uo 


|AiP 

Ui 




1 


4 

i=l 


( 10 ) 


We now minimize the thermodynamic potential ft^^ with 
respect to the two pairing fields Aq and Ai at half fill¬ 
ing and zero temperature T = 0. A global minimum at 
Aq = Ai = 0 corresponds to a phase without Cooper 
pairs, which is the normal phase, whereas a global mini¬ 
mum of at nonzero values for the pairing fields corre¬ 
sponds to a superfluid phase. We take the interactions at 
the two sublattices to be equal, Uq = Ui = U and map 
out the phase diagram as a function of the interaction 
strength U/t and chemical potential difference h/t, see 
Fig§ a). Without imbalance and interactions the sys¬ 
tem is in the normal state (fts). For a large enough in¬ 
teraction, but still without a population imbalance, the 
thermodynamic potential is minimized by nonzero and 
equal pairing fields Aq = Ai 7 ^ 0 , which we refer to as 
the SFo phase ( 1 ^ 2 )- This in contrast to the so-called 
TT-phase, which we call here the SF^^ phase, where Aq 
and Ai have opposite sign. For large enough imbalance 
h the thermodynamic potential is indeed minimized by 
nonzero Aq and Ai having opposite signs and the ground 
state of the system is the SF^^ phase. We find both an 
SFt^ phase where the pairing fields have equal magnitude 
((74), |Ao| = I All, and an SF^^ phase with unequal pair¬ 
ing fields (r^s). For even larger imbalances h the system 
enters the normal state again We find that most of 

the above phase transitions take place with the order pa¬ 
rameters changing discontinuously, suggesting first order 
phase transitions. The exceptions are between ft 2 and 
fts ai U = 2 t, and between (^4 and (^5 at [/ = 4t, where 
the pairing fields change continuously. In the above men¬ 
tioned SFo phase the two pairing fields Aq and Ai take 
the same value, which means a constant total pairing field 
through out the lattice. This corresponds to a homo¬ 
geneous superfluid phase, like in the usual BCS theory. 
Interestingly, in the SF^^ phases the pairing fields take 
different values on the different sublattices and the cor¬ 
responding phase is not a homogeneous superfluid phase. 

Because we neglected the momentum dependencies in 
the Hamiltonian in Eq.(|^ and consequently ended up 
with a thermodynamic potential without momentum in¬ 
tegrals it is possible to even find analytic expressions for 
the minima of the thermodynamic potential ft and the 
pairing fields, Aq and Ai, for which these minima are 



FIG. 3: (a) Zero temperature phase diagram. Dashed lines 

denote continuous phase transitions, all other transitions are 
first order, (b) Local minima of ft^^ as functions of h at 
different values of U. The different curves represent — Qe 
in Tab. [^ where the color coding is the same as in a). The 
dashed, dotted, and solid curves are for the cases of SFo, SFtt, 
and normal phases, respectively. 


acquired. For half filling, we list all possible local min¬ 
ima and the phase they correspond to in Tab.|T] together 
with the pairing fields. In Fig. |^b), we demonstrate 
the evolution of these local minima with the chemical 
potential difference h for different interaction strengths 
U. From this comparison the global minimum can be 
identified, which is how the phase diagram in Fig. [^a) 
was obtained. Also listed in Tab jl] for each phase are the 
polarizations P, where the polarization is the difference 
in densities between the two spin species divided by the 
total density, P = {n^ — The spin compo¬ 

nent densities can be calculated from the thermodynamic 
potential, Ua- = —dft/djicj- 
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H 2 

H 3 

H 4 

Q 5 

He 

SFtt (metastable) 

SFo {P = 0) 

N (P = 0) 

SF 

. (P = 1/2) 

N (P = 1) 

0<h<U/2-t 

0<h<U/2 

0 < h < t 

\h-t\< U/4 

0<h<U/2 

0 < t < h 

0<t<U/2 

0<t<U/2 

0 <U/2<t 

0<U/4<t 

0<t<U/4 

0<U 

CM 

II 

0 

<1 

Ao = yt/V4 - 

Ao = 0 

II 

0 

<1 

Ao = U/4 + yt/Vl 6 - P 

Ao = 0 

> 

II 

1 

to 

Ai = yt/V4 - ^2 

Ai = 0 

1 

II 

< 

Ai = -[//4+ yt/vi 6 -p 

Ai = 0 

II 

1 

to 

n = -2P/U - U/2 

Q = -2t 

^ = -h-t-U/8 

n = -h- 2P/U - U/4 

II 

1 

to 


TABLE I: Possible local minima of the thermodynamic potential in Eq.(lO) and their corresponding conditions at half 
filling, /i = 0, for equal interactions Uo — Ui — U. The polarization P is also shown for each phase. 


A. Connection between the inhomogeneous SF^r 
and LO phase 

Thus, we find from the above simplified, one¬ 
dimensional model with alternating signs for the hopping 
parameter, corresponding to a lattice with alternating s 
and p orbital sites, that in the presence of an imbalance 
between the two spin components the ground state of 
the system can be formed by an inhomogeneous super¬ 
fluid phase, the phase. Namely, in contrast to a 
homogeneous superfluid phase, the SFq phase, where the 
pairing field is constant throughout the system, in the 
SFt^ phase the pairing field changes with position in the 
lattice. 

We can compare the SF^^ phase to another inhomoge¬ 
neous superfluid phase, the so-called Larkin-Ovchinnikov 
(LO) phase [36]. In an LO superfluid the order parame¬ 
ter is also taken to be position dependent and specifically 
to be a cosine, where the wave vector q is left as a free 
parameter 

A(x) = Alo cos(gx). 

This Cooper pair ansatz results in a thermodynamic po¬ 
tential that depends both on the pairing field amplitude 
Alo and the wave vector g, (^(IAloI?^) [37]. Physically, 
the above Cooper pair ansatz corresponds to pairs formed 
by two fermions with different momenta, e.g., and 
The order parameter wave vector is equal to the 
net momentum q of the pairs and its wavelength is thus 
inversely proportional to it, Apo = 27r/g. 

Now, if the LO wavelength is twice the lattice spacing, 
Alo = 2 d, the LO phase in a lattice strongly resembles 
the SFt^ phase we find, where the pairing fields Aq and 
Ai only differ in sign (f^ 4 ), see FigJ^ Namely, the LO or¬ 
der parameter then takes the same value Apo, but with 
opposite sign, on neighboring sites. The SF^^ phase where 
the pairing fields also have a different magnitude can 
be viewed as a combination of a constant and a standing 
wave order parameter. In both cases, the SF^^ phase cor¬ 
responds to LO Cooper pairs with a net momentum of 
q = 7r/d, such as pairing with i^T^jd-kp The reason 
we can find this LO-like superfluid phase, without tak¬ 
ing it into account explicitly is because in a lattice the 
above pair corresponds to two particles with the same 
lattice momentum, k and —k, since 7r/d is the size of the 
Brillouin zone for a lattice where a unit cell contains two 



FIG. 4: LO order parameter and the SFtt phase order param¬ 
eters. For Q 4 the axes origin is at zero, whereas for Q 5 it is 
at some nonzero value. 


sites. 

It is possible that if a full LO ansatz is taken into ac¬ 
count for this system, a standing wave with a different 
wavelength Apo is found to be the ground state of the sys¬ 
tem. However, the general statement remains true, that 
a spin imbalance can result in inhomogeneous superfluid 
phase. 

Next we proceed to study a two-dimensional lattice, 
where there are two different p orbitals on the B sublat¬ 
tice, and we include the momentum dependence of the 
dispersions. In that case, taking a full LO ansatz into ac¬ 
count would be more involved. Interestingly, we are still 
able to find inhomogeneous superfluid phases in a rather 
simple manner, via the possibility of the SFt^ phase. 


IV. EXPERIMENTALLY REALIZABLE 
SYSTEMS WITH HYBRIDIZED ORBITALS 

In order to demonstrate how the phases revealed in 
the simple model can be observed experimentally, here 
we introduce a concrete two-dimensional lattice which 
has been realized recently by the group of Hemmerich 
[7] for ultracold bosonic atoms. In their experiment, a 
checkerboard lattice is created by two sets of orthogonal 
laser beams with shallow sites on one sublattice (A) and 
deeper sites on the other sublattice {B). By proper tun¬ 
ing of the relative depths of the two sublattices a system 
can be created where the lowest energy level of the A 
sites (the s band) are in resonance with the first energy 
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level of the B sites (the Px and Py bands). On all sites of 
the B sublattice the s bands are fully occupied, such that 
particles in the Px or py band can not relax to the lowest 
band. We study the possibility of superfluid phases for 
a two-component Fermi gas with a population imbalance 
loaded into such a lattice. To this end, we use a Hamilto¬ 
nian that includes hopping between nearest neighboring 
sites and attractive on-site interactions. We start by con¬ 
sidering the full lattice potential and arrive at an effective 
Hamiltonian by using a tight binding approximation. 


A. Tight Binding Approximation 

The lattice potential used in experiment can be de¬ 
scribed by 

V{x,y) = -Vo\cos{kox) + e*® cos(fco2/)Pi (H) 

where Vq is the average potential depth and ko is the wave 
vector determining the lattice spacing, ko = Trjd. In the 
following we take the recoil energy Er = li?ko/ {2nn) as 
the unit of energy, and the distance between two adjacent 
minima of A and B sites d as the unit of length. 

If the lattice potential Vq is large enough, a tight bind¬ 
ing model can well describe the properties of the sys¬ 
tem. In this approximation, the traps at shallower A 
and deeper B sites can be taken as harmonic potentials 
by expanding the lattice potential around each site as 

« -2Vb(l + cos6>) + V^klO + cos6>)(a;^ + y^) 

= Eo + ^nuj\{x‘^ + y"^), 

^ —2Vb(l — cos6>) + Vb/co(l — cosO){x‘^ -h y‘^) 

= Eo + ^}ko%{x‘^ + y'^), 

where the energy levels of the oscillators are E^^ = 
Eq'^ + kjujA.B (n +1/2). The degeneracy of the harmonic 
oscillator energy levels in two dimensions is n +1 with cor¬ 
responding parity (—1)’^. Based on these energy levels, 
the 5-band at the shallow lattice sites E^ is in resonance 
with the Px and Py bands on the deeper lattice sites E^ 
when EoEkujj^/2 = F^^ + 3/ico’^/2, which gives a relation 
between the lattice depth Vq and the phase 0 in Eq.(pT]). 
By numerically calculating the energy bands from the 
full lattice potential, we find for a lattice depth Vb = 10 
that the energy bands are in resonance for 0 ~ O.SSGtt, 
while from the harmonic oscillator energy levels one finds 
0 ^ O.SGOtt. This difference is small, which ensures that 
the tight binding is a good approximation. 

In the following, we focus on the above mentioned three 
orbitals, s on the A sublattice and p^ and p^ on the 


B sublattice in resonance with each other. The chemi¬ 
cal potentials of the system are chosen such that at low 
enough temperatures other bands are either fully occu¬ 
pied or empty and therefore play no role in our present 
study on superfluidity. 



X k 


FIG. 5: (a) Lattice potential in Eq. (11) with Vb = 10 and 
0 ^ 0.5567r in one unit cell, with the shallower A site at the 
origin and the deeper B sites at four corners. The bottom 
of the potential at B sites is —24.67, while it is —15.33 for 
A sites. The two arrows indicate nearest neighbor hopping, 
(b) Corresponding dispersions of single particle states along 

kx - ky 


The lattice potential in one unit cell and the corre¬ 
sponding band structure obtained by solving numeri¬ 
cally the single-particle Schrodinger equation are shown 
in Eig.[^ There, the lowest dispersion corresponds to the 
lowest energy band in the B sublattice. The first three 
bands above this band are the bands of interest, p^, 
and 5 , hybridized by the hopping. The even higher dis¬ 
persions correspond to even higher energy bands in the 
lattice. It can be seen that the three bands of interest 
are well separated from the other bands and considering 
their hybridization, the Hamiltonian term corre¬ 

sponding to the nearest-neighbour hopping can indeed 
be described by Eq. (§. 

If we also include next nearest neighbor hopping, we 
can write down a total single-particle Hamiltonian iLg in 
the basis 4']^ = 

= ( 12 ) 

k 


with the matrix 

















-K = 


(13) 


/e^ + cos kx cos ky 2it^^ sin kx 

I — 2 it^^ sin kx Cq + cos kx cos ky 
\ —2ity^ sin ky 2t^y sin kx sin ky 


2ity^ sin ky 
2t^y sin kx sin ky 
€q + 2tyy cos kx cos ky 


(a) to,. 


- 2-1012 

(b) 


- 2-10 1 2 

FIG. 6: Dispersions of the three hybridized states along (a) 
kx = ky and (b) along kx with = 0, for the lattice with 
Vb = 10 and 0 ~ O.SSGtt. The solid curves are exact dis¬ 
persions solved numerically, while the dashed curves are the 
htted dispersions of the reduced Hamiltonian with only two 
parameters eo and t. 


A B 

where the on-site energy offsets Cq ’ for A and B sites 
were added. 

The hopping coefficients and energy offsets can be ob¬ 
tained by fitting the dispersions obtained from diagonal¬ 
izing the above Hamiltonian to the exact dispersions cal¬ 
culated numerically. As an example, for the hopping pa¬ 
rameters we find ~ 0.0747 and for the energy 

offsets ~ —11.42 and Cq ~ —11.41, in the case of lat¬ 
tice parameters Vb = 10 and 0 ^ 0 . 5567 r, whereas the 
next nearest neighbor hoppings are at least three orders 
of magnitude smaller. Therefore, it is possible to use 
a reduced Hamiltonian with only nearest neighbor hop¬ 
pings t = and also, because are almost the 

same, they can be replaced by one parameter cq. If we 
now fit the dispersions obtained from the reduced Hamil¬ 
tonian with the exact dispersions, we find t ~ 0.0751 and 
eo ~ —11.42. The three dispersions of the new hybrid 

states are E^. = eo and = eo ± 

which reproduce the numerical results very well, see 
Fig-U In the following we will use the reduced Hamilto¬ 
nian. 

Here, it is worth pointing out that band E^ is exactly 


flat, i.e., dispersionless, which results from the linear 
combination of the and orbitals in the nearest- 
neighbor hopping approximation. Actually, the disper¬ 
sions depicted in Fig. are the same as the dispersions in 
the Lieb lattice, which illustrates the mapping discussed 
in the previous section. A final remark in this section 
is that the discussion above considers only the orbital 
degrees of freedom, and is valid for both fermionic spin 
species we consider. 

B. On-site Interactions 

To study the possibility of a pairing instability, like 
previously, we now include interactions and study the 
different superffuid phases that can occur in this two- 
dimensional system. 

We include an attractive s-wave contact interaction 

Hi = U j drip\{r)ipl{r)ipi{r)ip^{r), (14) 

where the interaction strength U < 0. By expanding the 
fermionic operators using the Wannier states the inter¬ 
action Hamiltonian reads, 

Y. j dr<(r-RX(r-R) 

= ^nin2n3n4'0R^/'0R^/'0R^^^ (1^) 

R,{nH 

with rii denoting the s, Px and Py orbitals and R the po¬ 
sition of the unit cell, and where we used the localizing 
property of the Wannier functions. The effective interac¬ 
tion coefficients I/m 712713714 absorb the corresponding cross 
integrals of four Wannier functions and are independent 
of R since all unit cells are equivalent in an infinite lat¬ 
tice. We use the harmonic oscillator eigenstates as an 
approximation to the Wannier functions to calculate the 
interaction coefficients. We only need to consider com¬ 
binations of the s band and the neighboring Px and Py 
bands within one unit cell, since all other cross integrals 
are at least four orders of magnitude smaller and can 
therefore be neglected. The results of the dominant in¬ 
teraction integrals are shown in Table |n| where it is used 
that the absolute value of a cross integral does not de¬ 
pend on the order of the Wannier functions 
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Uo — Ussss 

Ul - UxXXXf Uyyyy 

U2 — Uxxyy 

4.51 

4.04 

1.35 


TABLE II: Numerical values of Unin 2 n 3 n 4 ^lU for the lattice 
with Vo = 10 and 0 ~ O.SSGtt. 


The effective coupling constants are defined Uq = 


U.. 


Ui = u. 


= Uyyyy, aHcl U2 = U^cxyy, whcre 


in the harmonic approximation Ui = 8 ^/ 2 , indepen¬ 
dent of the lattice potential depth. In the U 2 interac¬ 
tion terms, the four orbitals yield six different combina¬ 
tions, namely , ip^^ipj^ip'^ip'^, 

and these terms with x and y interchanged, which are all 
included in our model. 

Using a simple mean-field approximation we introduce 
the BCS order parameters U,nnnm('0k;'0Ukt)- 

Considering the lattice symmetry, we have three differ¬ 
ent pairing fields denoted as Aq = A^^, Ai = A^^ = 
A^^, and A 2 = A^^ = A^^. For example, for the 
interaction term the mean-field approxi¬ 
mation is as follows 


- U2{^Hf){')Pl^) 




3Ui 


(16) 


For all other interaction terms the mean-field approxima¬ 
tion is similar. 


With the mean-field pairing, we understand that there 
are four interacting channels included in our Hamilto¬ 
nian, namely (n = s,x,y) counts intra¬ 

band pairing, ip^^ipj^^ip'^ip'^ yields interband pair tun¬ 
neling, ip^^ip^^'p'^ip^ results in interband pairing and 

ip^^ip^^ip^ip^ corresponds to spin exchange within inter¬ 
band pairs. However, as shown below, the last two inter¬ 
band pairing terms turn out to have no contribution. 


C. Full Hamiltonian 


Including the nearest-neighbour hopping and the pair¬ 
ing terms, as well as a population imbalance, the total 
mean-field Hamiltonian can be written with the Nambu 
basis 


H 

N 


■- {^k^BCS^k + 3[eo - (Ai - h )]} 


|Aop 

Uo 


8|Aip 

3Ui 


4|A2|2 
U 2 ’ 


(17) 


with the matrix 


^BCS = 


1 

1 

0 

2 it sin kx 

2 it sin ky 

Ao 

0 

0 

— 2 it sin kx 

1 

1 

0 

0 

0 

4Ai/3 

2 A 2 

—2it sin ky 

0 

1 

1 

0 

0 

2 A 2 

4Ai/3 

A* 

0 

0 

-eo + y - h 

— 2 it sin kx 

— 2 it sin k, 

0 

4A^/3 

2A5 

2it sin kx 

-eo + y-h 

0 

\ 0 

2AJ 

4Ap3 

2it sin ky 

0 

—eo y — 


(18) 


Accordingly, the thermodynamic potential reads 


^(^0, Ai, A 2 ) — — 


k 

1 


3[eo - (/i - h)] 


-V 7 > In 




1 + e 


-yuji(k) 




|Ao|2 8|Ai|2 4IA2I 


UqV 3 UiV U2V 


(19) 


where V is the 2D volume of a unit cell, C(;^(k) are the six 
eigenvalues from the 6 x 6 matrix in Eq. ( [l^ , and the 
quasi-momentum summation is over the first Brillouin 
zone. 


RESULTS 


From the thermodynamic potential Eq.(19) we can. 


like earlier, obtain phase diagrams for the 2D lattice with 
s and p orbital sites. We obtain phase diagrams as func¬ 
tion of interaction U and imbalance h for different tem¬ 
peratures and lattice parameters by minimizing D with 
respect to the order parameters Aq, Ai and A 2 . For the 
parameter regimes we considered, we find that A 2 is al¬ 
ways zero by numerically minimizing D. Subsequently, 
we calculate the momentum distributions for the spin- 
particles for the different phases we find. But first we 
take a look at the dispersions. 
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FIG. 7: Dispersions Ui for (a) the SFo phase and for (b) 
the SFtt phase as functions of kx at hxed ky = 0, with the 
parameters t ^ 0.075 and eo = ^ —11.42, and lattice pa¬ 

rameters Vo = 10 and 0 ~ O.SSGtt. In panel (a) Aq ~ 0.187 
and Ai 0.185, which minimize D at h = 0, 1/ = 0.10 and 
T = 0.01, while in panel (b) the sign of Ai is reversed to make 
a comparison between the SFo and SFtt phases. 


A. Dispersions 

The thermodynamic potential is calculated from the 
eigenvalues of Mbcs in Eq.(p^, which in the nor¬ 

mal state (Aq = Ai = 0) correspond to the particle 
dispersions and in the case of a superfiuid phase to the 
quasi-particle dispersions. As in the ID case, we find su¬ 
perfiuid phases with both Aq and Ai being nonzero, hav¬ 
ing either the same sign (SFq) or the opposite sign (SF^^). 
In Figj^ the dispersions uji for the SFq phase, panel (a), 
and the SF^^ phase, panel (b), can be compared, where 
all parameters were chosen the same to calculate these 
figures and the only difference is an added minus sign to 
Ai for FigJ^b). The flat dispersions are the same for the 
two cases. 

The dispersions are shown for the population balanced 
system. In the presence of a population imbalance, i.e. 

> 0, the dispersions are shifted downwards (upwards) 
for the spin up(down) particles, which are then the ma- 
jority(minority) particles. Intuitively, it can then be un¬ 
derstood from these dispersions that depending on the 
imbalance it is energetically more favourable to either 
occupy quasi-particle states corresponding to the SFq 
phase or to the SFt^ phase. However, to obtain the exact 
phase diagram of course the full thermodynamic poten¬ 
tial should be minimized, which is what we do next. 


B. Phase diagrams 


behaviour as we found for the simple ID model. 

In Figj^ phase diagrams are shown for the lattice pa¬ 
rameters Vo = 10 and 0 ^ 0.5567r at different temper¬ 
atures. White regions correspond to the normal phase, 
red to the SFq phase and (darker and lighter) blue cor¬ 
responds to SFt^ phases. In contrast to the ID case, here 
the SFt^ phase with Aq = —Ai (fl 4 ) is missing or at least 
highly reduced. The SF^^ phase with |Ao| ^ |Ai| (fls) 
splits into two phases, one with |Ao| < |Ai| (SF^) and 
one with |Ao| > |Ai| (SF^). The split of this phase 
was to be expected, since now Uq is not equal to Ui and 
thereby the degeneracy between the two local minima of 
O is lifted. 

We also observe that, with increasing temperature, the 
superfluid phases shrink towards the larger U and smaller 
h corner, with the SF^^ phase completely disappearing for 
high enough temperatures. 

Furthermore, to study the effect of different hopping 
parameters t, we change the lattice potential via Vq and 
6>, which can also be modified experimentally. In this way, 
the parameters t and eo obtained from the fitting in the 
tight binding model, as well as the effective interactions 
Ui^ are modified. We consider both a shallower lattice 
with Vb = 8 and a deeper lattice with Vb = 12, where the 
phase diagrams for these two cases are plotted in Fig. Ha) 
and (b) respectively. The lattice with Vb = 8 corresponds 
to a larger hopping coefficient t ^ 0.110 than previously, 
meaning that the other energy scales in the Hamiltonian, 
the interaction U and the imbalance h, become effectively 
smaller. The result is that the same phases as before 
now occur for larger U and /i, which can be observed 
in the phase diagram Fig. [^a). The deeper lattice with 
Vo = 12 corresponds to a smaller hopping t ~ 0.0517 and 
we observe the opposite effect. The superfluid SF^^ phase 
region now shifts towards smaller h and U. 


Momentum distributions 


As a possible experimental signature of the SFq and 
SFt^ phases, we present the quasi-momentum distribu¬ 
tions of the particles which can be observed experimen¬ 
tally. Since in such experiment, the original spin particles 
rather than the quasi-particles are observed, the particle 
distributions should be obtained by rotating the quasi¬ 
particle basis back to the original particle basis. This 
can be carried out by using the eigenvectors of Mbcs in 
Eq. (18), which form the transformation matrix S that 


diagonalizes Mbcs^ Then, the particle occupation num¬ 
ber of the state reads 


We now present phase diagrams as functions of chem¬ 
ical potential difference h and interaction strength U. 
Here, U is the full interaction strength from which the 
effective interactions Uq and Ui are calculated and is dif¬ 
ferent from the interaction coefficient used in the ID case. 
By minimizing O with fixed /i = eo, we find numerically 
that A 2 always vanishes, while Aq and Ai have similar 


3 


(k) ’ 


( 20 ) 


where cjj(k) are the eigenvalues of in Eqj^ 

Within the six bands given in the above equation, there 
are only three independent distributions, since + 
nk 4 , = 1 at half filling. Besides, since it is not possible to 
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FIG. 8: Phase diagrams as functions of h and U for lattice parameters Vb = 10 and 0 ~ O.SSGtt at different temperatures. In 
all phase diagrams the white region denotes the normal phase, red corresponds to the SFo and blue to the SFtt phases. The 
crosses mark the values for which the momentum distributions are calculated in Fig|10|and Fig|ll[ 



FIG. 9: Phase diagrams as functions of h and U for fixed tem¬ 
perature T = 0.01 at different lattice parameters. A shallower 
lattice is considered in (a) with Vb = 8 , ^ ~ 0.5647r, t ~ 0.110, 
eo ~ —8.44, Uq ^ 3.9817 and Ui ^ 3.6bU. A deeper lattice is 
used in (b) with Vb = 12 and 0 ~ O.SSOtt, with t ~ 0.0517, 
eo ~ —14.5, Uo ~ bU and Ui ~ 4.391/. As before, the white 
regions denote the normal phase (N), red corresponds to the 
SFo and blue to the SFtt phases. 


distinguish the original s and p bands of each particle in 
a measurement, we present the averaged occupations of 
the three bands for spin-up particles and for spin-down 
particles. 

Fig (T q| shows the quasi-momentum distributions for the 
points marked in FigJ^a), with interaction U = 0.10, 
temperature T = 0.01 and various values of the chemical 
potential difference h. The distributions show smooth 
changes as a function of the momenta. At zero tempera¬ 
ture, these would be sudden jumps corresponding to the 
Fermi surfaces of the filled bands. It can be seen that 
the momentum distributions are very different for the 
different phases. Especially, the difference between the 
SFo and SF^^ phase is considerable. Also the effect of 
the chemical potential difference on the densities can be 
seen quite clearly. At low /i, in Figjl^a), the differences 
between the t and ^ distributions are very small, mean¬ 
ing that the densities are similar. In contrast, at a large 
value of /i, in Figp^d), the difference between the t and 
^ distributions is very large, corresponding to a large po¬ 
larization. The various shapes of the distributions result 
from the interplay of the dispersions, being different for 
the various phases (Fi g0, the occupations of those lev¬ 
els, which depend on the chemical potential difference /i, 
and the temperature. 

The momentum distributions for the same interaction 
U = 0.10, but at a higher temperature T = 0.03 are 
shown in FigITTI for various values of the chemical po¬ 
tential difference /i, marked in FigJ^b). It can be seen 
that the qualitative differences between the momentum 
distributions for the different phases are still there, al¬ 
though a bit smoothened out compared to the T = 0.01 
distributions. However, the variations in the momentum 
distributions are now much larger, making the experi¬ 
mental observation of these interesting phases possible. 

At even higher temperatures the variations in the dis¬ 
tributions are even larger, but the qualitative differences 
between them for the various phases are then completely 
smoothened out. The quasi-momentum distributions 
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FIG. 10: Momentum distributions of the spin-up (top) and spin-down (bottom) particles averaged over the three bands as 
functions of kx and ky in a lattice with parameters Vo = 10 and 0 ^ O.SSGtt for the points marked in Fig ^a). The temperature 
is T = 0.01, the interaction is 1/ = 0.10 and the chemical potential differences (a) h = 0.10 (SFq), (1^ h = 0.14 (SF^), (c) 
h = 0.17 (SF^) and (d) h = 0.22 (normal state). 






but at temperature T = 0.03, corresponding to the points in Figl^b). Again the interaction 
I potential differences are (a) h = 0.10 (SFo), (b) h = 0.135 (SF^), (c) h = 0.17 (SF^) and (d) 

h — 0.22 (normal state). 
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FIG. 11: The same as in Fig. 
is U = 0.10 and now the chemica] 


also change with varying the interaction strength U. For 
the same phases the distributions are qualitatively the 
same as the ones depicted in Fig(TQ|and in FigfTT] How¬ 
ever, the variations in the distributions for both spin 
components become smaller at larger interaction strength 
and larger for smaller interactions U. 


VI. CONCLUSION AND OUTLOOK 

In conclusion, we studied lattices populated by two- 
component fermions occupying both s and p orbital 
states in both one and two dimensions. We showed how 
the system in two dimensions can be mapped to a Lieb 
lattice. In ID we used a simple mean-field calculation 


without including the full dispersions of the particles and 
determined the phase diagram, which shows two differ¬ 
ent superfluid phases. One superfluid phase is a homo¬ 
geneous superfluid phase, while the other one is an in¬ 
homogeneous superfluid phase, so-called tt phase, having 
similarities with the LO superfluid phase. Consequently, 
we calculated the full thermodynamic potential for an 
experimantally realizable two-dimensional lattice within 
a mean-field theory and find a similarly rich phase dia¬ 
gram. Also, we calculated the momentum distributions 
for the two spin components in the system, which could 
be observed experimentally. 

Due to the hybridization of the s and p bands a fiat 
band appears in the system. Flat bands can be related to 
many topological properties [381440] and may be respon- 
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sible for high Tc surface superconductor [3]. In future 
research, we will focus on the flat dispersion entering in 
this theory and the role of a flat band on pairing insta¬ 
bilities. 
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